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Abstract 

A  two-dimensional  axisymmetric  model  of  a  tubular  solid  oxide  fuel  cell  (SOFC)  is  presented.  The  model  analyzes  electrochemistry  and  heat 
transfer  inside  the  cell  with  new  approaches  to  thermal  and  chemical  simulations  being  introduced.  For  the  reforming  process  a  new  approach  based 
on  traditional  hypotheses,  equilibrium  and  kinetic  theory,  is  proposed.  Thus,  there  is  not  a  single  assumption  for  the  reaction  but  both  situations 
are  taken  as  possible  depending  on  composition  and  temperature;  results  show  that  equilibrium  is  reached  near  the  entry  zone  of  the  cell  and 
the  reaction  is  kinetically  controlled  later.  Thermally,  a  rigorous  study  of  radiative  heat  exchange  between  the  anodic  gas  and  the  anode  surface 
is  included  in  order  to  evaluate  its  effect  on  global  performance.  Results  show  that  this  effect  can  be  neglected  in  most  cases  but,  under  certain 
operating  conditions,  deviations  of  up  to  1%  may  appear. 

Additionally,  other  improvements  have  been  done  with  respect  to  previous  works:  viscosity  and  electrical  conductivity  of  gases  depend  on 
concentrations;  evaluation  of  concentration  losses  includes  a  new  description  of  binary  diffusion  coefficients;  calculation  of  convective  heat 
transfer  coefficients  includes  thermal  entry  length  and  laminar/turbulent  flow;  heat  generation  in  the  inner  part  of  the  cell  accounts  for  the  heat 
release  due  to  Joule  effect. 

©  2006  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

The  electric  power  generation  system  will  undergo  severe 
changes  in  the  following  decades.  This  transformation  is  aimed 
at  improving  the  performance  of  power  generators  in  three  dif¬ 
ferent  ways:  (i)  more  efficiency,  (ii)  lower  emissions  and  (iii) 
proximity  to  final  consumer. 

Starting  from  the  latest  goal,  proximity  to  consumer,  it  must 
be  said  that  the  current  layout  of  the  electric  power  transmission 
grid  is  highly  inefficient  in  two  ways.  Nowadays,  long  distance 
transmission  lines  are  used  to  connect  power  plants  with  load 
centres,  what  inevitable  leads  to  electric  losses  as  high  as  7% 
of  transmitted  power.  In  addition,  these  transmission  lines  are 
coupled  to  high  capacity  power  plants.  When  one  of  these  plants 
is  put  out  of  service  for  some  reason,  generation  and  load  are 
unbalanced  and  the  situation  may  not  be  easy  to  sort  out.  For 
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all  these  reasons,  there  is  agreement  between  all  partners  of  the 
electricity  market  to  evolve  at  least  partially  towards  a  so  called 
“distributed  generation  system”,  characterized  by  the  integra¬ 
tion  of  generators  into  consumers’  environment.  Thus,  a  new 
scenario  is  foreseen  where  rated  power  and  emissions  of  gener¬ 
ators  lower  while  efficiency  increases. 

Along  the  twentieth  century,  maximum  efficiency  values 
between  50%  and  60%  (LHV  based)  have  been  achieved 
by  fossil-fuel  based  power  plants,  with  outputs  ranging  from 
50  MW,  diesel,  to  800  MW,  natural  gas  combined  cycle.  In  spite 
of  this,  the  efficiency  of  this  kind  of  generators  when  applied  to 
a  distributed  system  is  significantly  lower  than  50%  due  to  the 
reduction  in  rated  power.  A  new  technology  is  needed  and  fuel 
cells  are  among  the  most  promising  options. 

Fuel  cells,  particularly  solid  oxide  fuel  cells,  are  free  of  ozone 
destructive  emissions  as  there  is  in  principle  no  carbon  monoxide 
or  nitrogen  oxide  in  the  exhaust  gas.  In  fact,  only  carbon  dioxide 
and  water  steam  are  emitted,  while  traces  (i.e.  less  than  1  ppm) 
of  NO*  may  be  generated  in  residual  combustion  processes. 
These  two  species  are  greenhouse  gases  but  the  total  amount 
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Nomenclature 

D  diameter  (m) 

De ff  diffusion  effective  coefficient  (m2  s-1) 

E  nernst  potential  (V) 

F  Faraday’s  constant  (C  mol-1) 

hcy  convective  heat  transfer  coefficient  (W  m-2  K- 1 ) 

Ahf  enthalpy  of  formation  (J  mol-1) 

i  current  density  (Am-2) 

i{)  exchange  current  density  (A  m-2) 

kcd  conductive  heat  transfer  coefficient  (W  m- 1  K- 1 ) 

Kq q  equilibrium  constant 

n  molar  flow  (mol  s - 1 ) 

nQ  number  of  exchanged  electrons 

Nu  Nusselt  number 

p  pressure  (Pa) 

Pr  Prandtl  number 

q  heat  flow  per  unit  area  (W  m-2) 

Q  total  heat  flow  (W) 

r  rate  of  reaction  (mol  m-2  s  - 1 ) 

R  gas  constant  (J  mol- 1  K- 1 ) 

Re  Reynolds  number 

S  cell  active  area  (m2) 

STCR  steam  to  carbon  ratio 

t  thickness  (m) 

T  temperature  (K) 

Uf  fuel  utilization  factor 

v  axial  coordinate  (m) 

xte.  thermal  entry  length  (m) 

Av  slice  length  (m) 

[X]  concentration  of  X 

Greek  letters 

a  transfer  coefficient  (activation  loss) 

tfgas  gas  absortivity 

y  pre-exponential  factor 

s  emissivity 

p  resistivity  (Q  m) 

a  Steffan-Boltzmann’s  contant 

Subscripts 
a  anode 

act  activation 

c  cathode 

cd  conduction 

con  concentration 

cv  convection 

e  electrolyte 

eq  equilibrium 

f  formation 

h  hydraulic 

ohm  ohmic 

ref  methane  reforming  reaction 

shift  shift  reaction 

0  standard 


Superscripts 
0  bulk  flow 

*  active  site 


Table  1 

Reference  geometry 


Length  (m)  1.5 

Anode  outer  diameter  (mm)  22 

Anode  thickness  (pm)  100 

Electrolyte  thickness  (pm)  40 

Cathode  thickness  (mm)  2.2 

Metallic  interconnection  thickness  (pm)  85 

Active  area  (cm2)  834 


is  far  less  than  with  conventional  plants  due  to  the  increase  in 
efficiency. 

Fuel  cells  are  thought  as  excellent  complements  to  heavy 
duty  power  plants  in  the  mid  and  long  term  but  their  technology 
must  be  improved  in  many  ways  so  that  both  performance  and 
reliability  increase  while  cost  remarkably  decreases.  Control  of 
temperature  inside  high  temperature  fuel  cells,  as  well  as  fore¬ 
cast  of  mechanical  failures  caused  by  temperature  oscillations 
must  be  studied  in  detail.  To  do  so,  a  model  of  performance  is 
needed.  This  model  must  give  a  precise  description  of  every  elec¬ 
trochemical  and  heat  transfer  process  taking  place  inside  the  cell. 


4  air 


Fig.  1.  Schematic  view  of  power  plant. 
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Finally,  a  fuel  conditioning  system  is  needed.  Although  internal 
reforming  of  natural  gas  is  possible  when  dealing  with  high  tem¬ 
perature  fuel  cells,  part  of  the  fuel  flow  must  be  reformed  before 
it  enters  the  anode  of  the  fuel  cell. 

The  work  presented  here  covers  all  the  subjects  mentioned 
above  and  gives  a  very  valuable  tool  to  further  develop  the  tubular 
solid  oxide  fuel  cell  technology.  Although  the  model  is  fully 
flexible  and  applicable  to  tubular  cells  with  any  dimensions,  a 
hundred  watt  Siemens- Westinghouse  cell  has  been  considered 
as  base  case;  see  Table  1  [1].  In  addition,  a  prereformer  fed  with 
natural  gas  and  anode  exhaust  gas  is  included.  A  schematic  view 
of  the  plant  can  be  found  in  Fig.  1. 

2.  Fuel  conditioning 

2. 1 .  Production  of  hydrogen 

Natural  gas,  consisting  mainly  of  methane,  is  the  most  inter¬ 
esting  source  of  hydrogen  for  stationary  applications  of  SOFCs 
[2].  Methane  conversion  to  hydrogen  and  by-products  can  be 
done  through  different  techniques,  steam  reforming  being  the 
most  suitable  if  fuel  cells  are  involved.  In  this  case,  when  a 
stream  of  methane  mixes  with  another  one  with  a  high  content 
of  water  steam,  the  following  reaction  may  take  place: 

CH4  +  H20  ->  3H2+CO  (1) 

Eq.  (1)  describes  the  reforming  process  and,  since  it  is  endother¬ 
mic  (A/zf  =  206kJkmol_1),  is  usually  coupled  to  high  temper¬ 
atures  or  to  heat  addition.  Carbon  monoxide,  which  is  a  by¬ 
product  of  the  reforming  reaction,  is  further  oxidized  to  carbon 
dioxide  if  water  is  still  present  (Eq.  (2)) 

CO  +  H20  ->  H2  +  C02  (2) 

This  second  reaction  (Eq.  (2)),  named  water  gas  shift  reaction, 
is  slightly  exothermic  (A/if  =  — 41  kJkmol-1)  and  favours  the 
formation  of  hydrogen. 

Eq.  (1)  and  (2)  take  place  inside  the  anodic  duct  of  an  SOFC 
due  to  the  high  operating  temperature  of  the  cell  and  the  high 
concentration  of  water  steam  in  the  anodic  gas,  which  can  be  as 
high  as  50%.  This  steam  is  formed  when  hydrogen  is  oxidized 
according  to  the  following  exothermic  (A/jf  =  —242  kJ  kmol-1) 
reaction: 

H2  +  02  ->  H20  (3) 

Nevertheless,  not  only  Eqs.  (1)— (3)  are  possible  and  some  other 
undesirable  reactions  can  be  found.  In  particular,  decomposi¬ 
tions  of  methane  or  carbon  monoxide  may  take  place  through 
the  following  reactions: 


ch4- 

►  C  +  2H2 

(4) 

2CO  - 

-*  c  +  co2 

(5) 

which  are  endothermic  (Ahf  =  75  kJkmol-1)  and  exothermic 
(A/jf=  — 172. 6  kJ  kmol- 1 ),  respectively.  If  any  of  these  two  reac¬ 
tions  took  place  at  the  anode,  deposition  of  carbon  particles  over 
the  anode  surface  would  block  the  catalyst  and  dramatically 


Temperature  (K) 

Fig.  2.  Gibbs  free  energy  change  for  reactions  (1),  (2),  (4)  and  (5). 

reduce  its  activity;  in  such  a  case,  the  electrochemical  perfor¬ 
mance  of  the  cell  would  drop  down. 

Temperatures  around  950  °C  enhance  those  reactions  whose 
behaviour  is  more  endothermic,  i.e.  Eq.  (1),  while  decreasing  the 
rate  of  exothermic  reactions,  specially  Eq.  (4)  and,  to  a  lesser 
extent,  Eq.  (2).  Fig.  2  shows  the  variation  of  free  energy  for 
reactions  Eqs.  (1)— (2),  (4)— (5)  with  temperature  [3]. 

The  changes  in  Gibbs  free  energy  are  related  with  the  way  in 
which  a  certain  reaction  occurs.  In  other  words,  a  reaction  whose 
Gibbs  free  energy  change  is  positive  at  a  certain  temperature 
will  shift  towards  reactants.  Oppositely,  the  same  reaction  would 
shift  towards  products  if  the  change  in  the  Gibbs  free  energy  of 
the  system  were  negative  at  that  temperature.  The  latter  case  is 
known  as  a  spontaneous  process.  According  to  the  relationship 
between  Gibbs  free  energy  and  spontaneity,  Fig.  2  shows  that 
below  800  K,  neither  methane  reforming  nor  decomposition  are 
spontaneous  processes;  at  this  temperature,  carbon  monoxide  is 
decomposed  into  particles  and  carbon  dioxide  rather  than  oxi¬ 
dized  with  water,  Eqs.  (5)  and  (2),  respectively.  Despite  this, 
if  temperature  rises  above  1 100  K  (850  °C),  methane  reforming 
and  decomposition  become  spontaneous.  Eq.  (3)  is  spontaneous 
at  any  temperature. 

Finally,  the  concentration  of  steam  is  increased  in  order  to 
enhance  the  conversion  of  methane  into  hydrogen  and  not  into 
atomic  carbon.  By  doing  this,  and  according  to  Le  Chatelier’s 
principle,  the  equilibrium  of  the  reforming  reaction  is  shifted  to 
the  products. 

The  exhaust  gas  of  the  anode  is  the  most  interesting  source 
of  both  heat  and  steam  for  the  reforming  process  inasmuch  as 
its  temperature  is  somewhere  near  1000  °C  and  it  contains  up  to 
50%  (v)  of  water  steam.  This  gas  is  recirculated  and  mixed  with 
the  fuel  before  this  enters  the  cell. 

2.2.  Prereforming 

It  has  been  stated  that  internal  reforming  of  methane  is  pos¬ 
sible  in  solid  oxide  fuel  cells.  However,  if  the  anode  is  fed  with 
pure  methane  temperature  drops  down  significantly  in  the  first 
part  of  the  cell  due  to  the  endothermy  of  the  reforming  reaction. 
Later  on,  when  methane  is  fully  reformed  to  hydrogen,  the  heat 
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Fig.  3.  Effect  of  prereforming  on  temperature  of  tube  at  0.6  V  and  80%  fuel 
utilization. 


released  in  Eq.  (3)  is  no  longer  consumed  by  Eq.  (1)  and  tem¬ 
perature  increases;  it  must  be  taken  into  account  that  the  total 
heat  released  by  Eq.  (3)  doubles  that  consumed  by  Eq.  (1)  [4]. 
Besides,  the  initial  temperature  drop  affects  the  electrochemical 
performance  of  the  cell  negatively  and  gives  place  to  high  ther¬ 
mal  stresses,  which  can  ruin  the  cell  due  to  the  combination  of 
a  significant  temperature  gradient  along  the  tube  and  different 
thermal  expansion  coefficients  of  anode,  cathode  and  electrolyte 
[5].  Consequently,  methane  must  be  partially  reformed  in  an 
external  reactor  [2,4] . 

Fig.  3  shows  that  the  evolution  of  tube  temperature  along 
the  cell  is  smoother  when  the  degree  of  prereforming  is  high. 
Nevertheless,  this  situation  also  leads  to  a  drop  in  electrochem¬ 
ical  performance  so  an  intermediate  working  point  must  be 
adopted. 

2.3.  P re refomer  model 


K 


eq,  shift 


[H2][C02] 

[CO][H20] 


log  K& q  =  aT4  +  bT 3  +  cT2  +  dT  +  e 


Recirculation  of  anode  exhaust  flow  is  done  in  order  to  avoid 
the  problems  mentioned  before,  the  amount  of  anode  exhaust 
flow  to  be  recirculated  being  related  to  the  molar  flow  of 
atomic  carbon  entering  the  reactor  through  the  “steam  to  car¬ 
bon  ratio”  or  STCR.  STCR  is  the  ratio  of  steam  molar  flow  to 
methane  and  carbon  monoxide  molar  flow  at  the  prereformer 
inlet 


STCR  = 


^h2o 


^ch4  +  ^co 


inlet  prereformer 


According  to  Peters  et  al.  [8],  for  STCR  higher  than  2/2.5  carbon 
deposition  does  not  take  place  inside  the  prereformer. 

The  prereformer  is  considered  to  be  adiabatic,  i.e.  no  heat  is 
transferred  to  or  from  the  environment.  Thus,  the  composition 
of  the  exhaust  gas  can  be  calculated  by  applying  simultaneously 
the  laws  of  conservation  of  mass  and  energy  and  the  laws  of 
chemical  equilibrium  (Eqs.  (6)-(8))  at  the  exit  of  the  reactor. 


3.  Fuel  cell  model 

Fig.  1  shows  reference  geometry  for  the  base  case.  Air  enters 
the  injector  tube  and  flows  downwards.  At  the  end  of  the  tube, 
the  flow  turns  back  to  the  upwards  direction  through  the  annu¬ 
lar  duct  between  the  injector  tube  and  the  inner  wall  of  the 
cathode.  Parallel  to  it,  a  fuel  stream  coming  from  the  prere¬ 
former  flows  upwards  through  the  anode  duct  at  the  outer  face 
of  the  cell  tube.  This  layout  can  be  thought  of  as  coflow  and  is 
divided  into  a  number  of  slices  in  the  axial  direction.  Each  one 
of  these  elements  is  divided  again  into  five  elements  in  the  radial 
direction: 


The  prerefomer  is  a  catalytic  bed  reactor  fed  by  streams  of 
natural  gas  and  anodic  exhaust  gas.  These  gases  react  according 
to  Eqs.  (1)  and  (2)  which  are  considered  to  reach  equilibrium  at 
the  exit  temperature.  The  composition  of  the  effluent  is  defined 
by  constants  Ke q?ref  and  7feq, shift  which  depend  on  temperature 
according  to  the  polynomial  in  Eq.  (8).  Coefficients  in  Eq.  (8)  are 
different  when  applied  to  KQ q?ref  and  #eq, shift  and  can  be  found 
in  Table  2  [3,6,7] 

_  [H2]3[CO] 

eq,ref  [CH4][H20]  (  ’ 


Table  2 

Coefficients  in  Eq.  (8)  [3,7] 


Reforming 

Shift 

a 

-2.6312  x  1CT11 

5.47  x  10“12 

b 

1.2406  x  10“7 

-2.5748  x  10“8 

c 

-2.2523  x  10“4 

4.6374  x  10“5 

d 

1.9503  x  10-' 

-3.9150  x  10“2 

e 

-66.1395 

13.2097 

1.  Air  injector  tube. 

2.  Air  injector  tube  wall. 

3.  Annular  duct  (cathode). 

4.  Cathode-electrolyte-anode  wall. 

5.  Anode  duct. 

The  division  made  is  standard  for  a  finite  volume  method 
application  to  tubular  cells  and  has  been  used,  with  few  changes, 
by  Campanari  [1]  and  Stiller  et  al.  [9]  before.  Fig.  4  shows  a  slice 
and  its  elements. 

Axial  symmetry  is  considered,  both  for  the  electrochemical 
and  thermal  models,  on  the  basis  of  the  distribution  of  tubes 
inside  a  fuel  cell  stack  as  shown  in  Fig.  5.  A  stack,  in  particular  for 
a  Siemens- Westinghouse  cell,  is  formed  by  a  number  of  bundles 
of  24  tubular  cells  on  a  3  x  8  configuration.  Those  cells  located  in 
the  core  of  the  stack,  A  in  Fig.  5,  are  surrounded  by  other  cells 
working  under  similar  conditions.  In  such  cases,  considering 
axial-symmetric  boundary  conditions  is  acceptable  and  the  polar 
coordinate  can  be  eliminated  from  the  calculations.  If  the  cell  is 
located  in  the  outer  part  of  the  stack,  B  in  Fig.  5,  axial- symmetric 
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i-1 


SLICE  i 


i+1 


4.  Internal  reforming 


Fig.  4.  Elements  in  a  slice. 

boundary  conditions  cannot  be  considered  and  a  fully  three- 
dimensional  analysis  must  be  done. 

Fuel  cell  stacks  are  usually  formed  by  a  lot  of  tubes,  say  1000 
tubes  per  kilowatt,  and  most  of  them  are  located  in  its  core  so  a 
hypothesis  of  axial  symmetry  is  applicable  to  most  of  the  tubes 
in  a  stack.  However,  if  the  effects  of  peripheral  tubes  are  to  be 
considered,  some  authors  consider  a  2-3%  penalization  in  total 
power  due  to  heat  losses  to  the  surroundings  [10]  according  to 
Siemens- Westinghouse  recommendations . 

Finally,  in  order  to  decrease  the  computational  load,  the 
analysis  of  the  cell  is  split  into  two  blocks.  Firstly,  the  cell 
electrochemical  behaviour  is  studied  and,  for  this  purpose,  a 
prescribed  temperature  condition  is  set  so  that  the  temperature 
field  is  invariable.  Later,  an  analysis  of  heat  transfer  phenomena 
inside  the  fuel  cell  is  done  where  current  density,  voltage  loss, 
reaction  rate  distributions  are  known.  Although  both  blocks  can 
be  solved  independently,  they  are  inseparables  from  a  thermody¬ 
namic  point  of  view.  Fig.  6  shows  that  the  solution  to  the  analysis 
is  single. 


Fig.  5.  Bundles  in  a  fuel  cell. 


Previous  to  describing  the  fuel  cell  thermal  and  electrochemi¬ 
cal  models  of  performance  some  comments  on  methane  internal 
reforming  must  be  done.  The  anodic  duct  of  a  solid  oxide  fuel 
cell  is  fed  by  the  effluent  of  the  prereformer  which  consists 
mainly  of  hydrogen  and  water  and  has  about  10-20%  (v)  of 
unreformed  methane  and  carbon  monoxide.  As  long  as  methane 
is  still  present  two  heats  of  reaction  tend  to  heat  up  and  cool  down 
the  circulating  gas:  heat  released  by  the  oxidation  of  hydrogen 
and  consumed  by  the  reforming  of  methane.  Carbon  deposi¬ 
tion  is  again  avoided  due  to  the  high  temperature  and  content  of 
steam. 

The  water  gas  shift  reaction  is  considered  to  reach  equilib¬ 
rium,  defined  by  constant  shift  (Eq.  (7)),  at  each  point  of  the 
duct.  The  temperature  to  be  considered  in  Eq.  (8)  is  that  of  the 
anodic  gas  at  each  point  of  the  duct. 

For  the  reforming  reactions  two  scenarios  are  considered 

1 .  Methane  reforming  reaches  equilibrium.  Under  these  circum¬ 
stances,  the  analysis  is  similar  to  that  of  the  prereformer. 
Eqs.  (6)— (8)  along  with  heat  and  mass  balances  are  used  to 
determine  the  equilibrium  concentration  at  the  anodic  local 
temperature. 

2.  Methane  reforming  reaction  is  slow  and  equilibrium  is  not 
reached.  The  process  is  kinetically  controlled.  Ahmed  and 
Fogger’s  [11]  equation  for  the  reforming  rate  is  used  in  this 
case  (Eq.  (10)) 

_  7  0,85  0,35  (  95  xl03^  /im 

rCH4  —  ^refPcH4PH20  eXP  I  \ 

It  must  be  pointed  out  that  both  hypotheses  have  been  used 
before  to  describe  the  reforming  process  inside  a  high  temper¬ 
ature  fuel  cell.  Thus,  Magistri  et  al.  [12]  assume  equilibrium 
conditions  while  Campanari  [1],  Stiller  et  al.  [9]  or  Selimovic 
[7]  employ  kinetic  theory  equations.  However,  the  methodology 
presented  in  this  work  is  a  new  approach  to  the  internal  reform¬ 
ing  analysis  as  long  as  it  takes  into  account  that  this  process  can 
be  characterized  by  different  mechanisms  in  different  positions 
inside  the  cell. 
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The  following  method  is  used  to  evaluate  the  reforming  rate 
of  reaction.  First,  the  rate  of  reaction  is  calculated  considering 
that  equilibrium  is  reached.  Then,  this  value  is  compared  to  that 
given  by  Eq.  (10)  after  an  iterative  calculation  to  include  the 
effect  of  water  gas  shifting.  If  the  rate  of  reaction  given  by  the 
kinetic  theory  is  higher  than  that  at  equilibrium,  equilibrium  is 
reached.  Else,  equilibrium  is  not  reached  and  Eq.  (10)  must  be 
used.  Fig.  7  shows  the  rates  of  reaction  when  each  of  the  models 
proposed  are  used:  equilibrium,  kinetics  and  both.  For  the  latter, 
the  “mixed  method”  denomination  has  been  used. 

Tree  different  zones  can  be  identified  in  Fig.  7.  The  first  zone 
is  marked  A  and  ranges  from  0  to  30  cm  in  the  axial  direction  of 
the  fuel  flow.  In  this  zone,  the  reforming  process  reaches  equi¬ 
librium  and  the  rate  of  reaction  is  determined  by  the  equilibrium 
constant  Keq  ref.  Oppositely,  in  zone  C  the  process  is  controlled 
by  kinetics  because  the  rate  of  reaction  is  too  slow  due  to  a  low 
concentration  of  methane  and  equilibrium  is  not  reached;  i.e.  the 
rate  of  reaction  assuming  equilibrium  is  far  higher  than  that  pre¬ 
dicted  by  Eq.  (10).  Zone  B  is  a  transition  zone  and  temperature 
determines  whether  equilibrium  is  reached  or  not. 

Fig.  7  shows  that  the  rate  of  reaction  when  the  “mixed 
method”  is  used  is  close,  slightly  above,  to  the  rate  predicted 
by  an  equilibrium  hypothesis  initially  and  close,  again  slightly 
above,  to  the  rate  of  reaction  calculated  from  Eq.  (10)  later  on. 
The  reason  why  the  “mixed  method”  line  does  not  coincide  with 
the  equilibrium  and  kinetics  lines  in  zones  A  and  C,  respectively 
is  that  fuel  and  air  flow  in  opposite  directions  so  temperatures 
at  one  end  of  the  tube  strongly  depend  on  temperatures  at  the 
other  end.  The  influence  of  the  reforming  description  on  tube 
temperature  is  shown  in  Fig.  8. 

Table  3  shows  the  effect  of  different  descriptions  of  the 
reforming  process  on  final  results.  Two  conclusions  can  be 
drawn.  The  maximum  error  made  when  using  any  of  the  two 
usual  descriptions,  kinetics  or  equilibrium,  is  1%  on  the  main 
parameters  of  performance.  Besides,  Table  3  shows  that  assum¬ 
ing  the  reforming  reaction  at  equilibrium  underestimates  power 
and  efficiency  while  using  Eq.  (10)  leads  to  an  overestimation, 
in  both  cases  with  respect  to  the  “mixed  method”,  which  is  con¬ 
sidered  to  be  closer  to  reality  as  shown  later. 


Table  3 

Effect  of  reforming  description  on  results  (pressure  =  3  bar) 


Kinetics 

Equilibrium 

Mixed 

Voltage  (V) 

0.6 

0.6 

0.6 

Current  density  (Am-2) 

3045 

2988 

3002 

Fuel  utilization  (%) 

81.2 

79.7 

80 

Air  utilization  (%) 

19.3 

18.9 

19 

Power  (W) 

152.4 

149.5 

150.2 

Efficiency  (%) 

46.44 

45.56 

45.79 

Mean  temperature  (°C) 

892 

895 

896 

Peak  temperature  (°C) 

1010 

999 

1003 

5.  Electrochemical  model 


According  to  the  discretization  shown  above  (Fig.  4),  an 
electrochemical  analysis  based  on  deviations  from  ideal  per¬ 
formance  is  proposed  at  each  slice.  Considering  the  electrodes 
to  be  equipotential  surfaces  [10,13],  the  ideal  or  Nernst  voltage 
is  defined  by 


E  =  Eo  + 


RT 


In 


/  1  : 
I  P^iPo2 

\  Pu2o 


(11) 


where  Eo  is  the  standard  potential,  which  varies  almost  lineally 
with  temperature  [3] 

Eq  =  1.26485  -  2.4725  x  10~4T  -  1.875  x  10 ~8T2  (12) 


Irreversibilities,  also  known  as  polarizations  or  overpotentials, 
arise  that  cause  voltage  losses  from  the  Nernst  voltage  value  and 
have  a  negative  impact  on  fuel  cell  performance.  There  are  three 
types  of  overpotentials: 


1 .  Activation  loss :  A  Vact . 

2.  Ohmic  loss:  AV0hm. 

3.  Concentration  loss:  AVcon. 


Multidimensional  models  previously  developed  frequently 
include  activation  and  ohmic  losses  but  not  concentration  loss 
as  this  is  one  order  of  magnitude  smaller  under  usual  operating 
conditions. 
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5.1.  Activation  overpotential 

Activation  losses  take  place  at  the  electrodes  and  are  related 
to  the  minimum  voltage  drop  needed  to  overcome  the  activation 
energy  of  both  hydrogen  and  oxygen  ion  forming  reactions.  The 
activation  overpotential  for  each  electrode  is  usually  defined  by 
the  Butler- Volmer  equation 


complexity  of  evaluating  charge  transfer  coefficients  is  applica¬ 
ble  to  the  coefficients  involved  in  Eqs.  (18)  and  (19). 

5.2.  Ohmic  overpotential 

Ohmic  overpotential  is  caused  by  the  flow  of  electric  charge 
across  electrodes  and  electrolyte.  It  can  be  defined  by 


i  =  i  o 


TIqFA  Vact  \ 
RT  ) 


—  exp 


n&FA  VaCt  \  1 

FT  Jj 


(13) 


where  i o,  aa  and  ac  are  the  exchange  current  density  and  transfer 
coefficients,  respectively. 

Charge  transfer  coefficients  are  related  to  the  relative 
response  of  current  density  at  each  electrode  to  a  change  in  acti¬ 
vation  overpotential.  Evaluating  these  coefficients  is  a  difficult 
task  and  they  are  usually  considered  symmetric  and  equal  to  0.5 
for  a  one-electron  charge- transfer  reaction  [14].  In  such  a  case, 
Eq.  (13)  can  be  rewritten  as 


i  =  2/o  sinh 


f  KeFAVact  \ 

V  ft  J 


(14) 


Eq.  (14)  can  be  further  simplified  assuming  that  activation  losses 
are  high,  Tafel  equation  (Eq.  (15)),  or  low,  linear  equation  (Eq. 
(16)) 


A  Vact  =  a  +  b  In  i 


RT 

A  Vact  =  —=ri 

n&Fi0 


(15) 

(16) 


Activation  losses  are  considered  low  if  Eq.  (17)  is  satisfied  [7] 


FAVact 

neRT 


(17) 


However,  this  symmetric  relationship  of  charge  transfer  coeffi¬ 
cients  is  not  based  on  any  physical  law  and  the  only  accurate  way 
of  determining  aa  and  ac  is  to  do  some  testing.  In  fact,  it  is  not 
even  clear  whether  transfer  coefficients  are  constant  for  a  given 
fuel  cell  and  reaction  or  they  depend  on  current  density.  For  these 
reasons,  data  from  tests  done  by  Costamagna  and  Honegger  [15] 
on  SOFCs  have  been  used  and  simplifications  of  Butler- Volmer 
equation  have  been  discarded  in  this  work.  It  must  be  said  that 
calculating  activation  overpotentials  from  simplified  equations 
can  lead  to  errors  as  high  as  23%  as  shown  by  Chan  et  al.  [16]. 

The  exchange  current  density  is  related  to  the  current  density 
exchanged  back  and  forth  when  equilibrium  is  reached  and  is 
not  easy  to  evaluate  either  [14].  Despite  this,  there  is  agreement 
that  the  following  reactions  are  appropriate  for  anode  (Eq.  (18)) 
and  cathode  (Eq.  (19)): 


*0,a  —  Va 


Ph2\ 
Pref  J 


)o 

exp 


F  act,  a 

RT 


^0,c  —  Yc 


)c 

exp 


F  act,c 

RT 


(18) 

(19) 


Pre-exponential  factors,  exponents  a,  b  and  c  and  activation  ener¬ 
gies  are  taken  from  tests  by  Costamagna  and  Honegger  [15].  The 


AVohm  =  iFQ  f  (20) 

where  Re f  is  taken  as  an  effective  electric  resistance  similar  to 
that  of  the  structure  formed  by  electrolyte,  electrodes  and  inter¬ 
connection  plates.  Complexity  of  current  flow  paths  inside  the 
cell  makes  it  difficult  to  evaluate  an  equivalent  resistance.  In 
this  work,  the  model  developed  by  Nisancioglu  [17],  and  pre¬ 
viously  used  by  Haynes  [10],  Campanari  [1]  and  Stiller  et  al. 
[9],  has  been  completed  with  a  temperature  dependence  of  elec¬ 
tronic/ionic  resistivities  which  was  not  included  in  [17].  This 
dependence  is  taken  from  ref.  [7]  and  can  be  expressed  by  the 
following  equation: 

T  ( B\ 

Pi(T )  =  -exp  (-)  (21) 

where  A  and  B  depend  on  materials. 

5.3.  Concentration  losses 


Concentration  overpotential  is  related  to  the  velocity  at  which 
diffusion  of  species  takes  place  inside  porous  electrodes.  Gen¬ 
erally  speaking,  it  can  be  said  that  there  are  two  main  targets 
for  electrode  designers.  First,  it  is  important  that  electrodes  be 
thin  because  the  thinner  the  layer  through  which  species  diffuse, 
the  faster  the  hydrogen  oxidation  reaction  and  the  higher  the 
current  density  of  the  cell.  The  lower  limit  of  electrode  thick¬ 
ness  is  determined  by  available  manufacturing  techniques  and 
mechanical  integrity  matters.  Besides,  operating  conditions  have 
a  notorious  impact  over  concentration  polarization.  In  particu¬ 
lar,  species  flows  at  each  electrode  are  caused  by  concentration 
gradients  between  bulk  flow  and  active  sites.  For  this  reason, 
when  one  of  the  reactants  is  near  to  exhaustion,  the  whole  cell 
performance  is  controlled  by  the  slower  diffusion  process  what 
causes  a  higher  voltage  drop.  This  situation  can  also  be  found  if 
any  of  the  reactants  is  near  to  100%  conversion  or  if  current  den¬ 
sity  becomes  too  high. According  to  Selimovic  [7],  concentration 
overpotential  can  be  evaluated  from  the  following  expression: 


AVcon 


*Li n  ( +  5Lln  ( Ek) 

n*F  \Ph2Pu2  o)  n*F  UoJ 


(22) 


where  superscripts  “0”  and  stands  for  properties  at  free  bulk 
flow  and  active  sites,  respectively.  If  concentration  losses  are 
to  be  known,  properties  at  active  sites  must  be  related  some¬ 
how  to  free  flow  properties  as  long  as  only  these  are  known 
when  analysing  the  cell.  Fick’s  laws  are  used  to  evaluate  con¬ 
centrations  at  active  sites  from  free  flow  concentrations.  For  this 
purpose,  effective  diffusion  coefficients  are  developed. 

The  calculation  of  diffusion  coefficients  is  out  of  the  scope  of 
this  paper  and  would  not  be  described.  Still,  it  must  be  said  that 
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calculating  diffusion  coefficients  of  gases  which  are  part  of  either 
cathodic  or  anodic  flows  is  a  difficult  task.  Experimental  data 
are  very  scarce  and,  where  available,  they  are  taken  at  very  low 
temperatures  compared  to  that  inside  an  SOFC.  This  situation 
gets  worse  when  applied  to  binary  diffusion  coefficients. 

A  solution  to  this  problem  has  usually  been  found  in  using 
semi-empirical,  theoretically  based  methods,  which  extend 
available  data  to  practical  operating  conditions.  Lennard  Jones’ 
method  for  binary  coefficients  is  the  most  popular  and  extended, 
although  Todd  and  Young  [18]  show  that  using  this  method  leads 
to  remarkable  deviations  from  experimental  data  at  high  temper¬ 
ature.  Instead,  Fuller  et  al.’s  method  is  shown  to  be  much  more 
accurate,  especially  at  high  temperature,  what  makes  it  suitable 
for  fuel  cell  modelling. 

Eq.  (23)  and  (24)  are  obtained  when  applying  Fick’s  laws  as 
described  in  [19]  and  represent  concentration  losses  at  cathode 
and  anode,  respectively 


Table  4 


Input  parameters  to  model 


Anode  resistivity  (£2  m)  (Eq.  (21)) 

pa  =  (779 5  X  106)  exp  (- 1 150/D 

Cathode  resistivity  (Q  m)  (Eq.  (21)) 

pc  =(7742  x  106)  exp  (- 1200/D 

Electrolyte  resistivity  (£2m)  (Eq.  (21)) 

pe  =  (1/3.34  x  104)exp(- 10300/D 

Interconnection  resistivity  (Q,  m)  (Eq. 
(21)) 

picm  =  (779.3  x  106)  exp  (- 1100/D 

Anode  conductivity  (W  m-1  K-1) 

3 

Cathode  conductivity  (W  m-1  K-1) 

3 

Electrolyte  conductivity  (W  m- 1  K- 1 ) 

3 

Anode  emissivity 

0.9 

Cathode  emissivity 

0.9 

Butler- Volmer,  anode  (Eq.  (13)) 

aa  =  2,  ac  =  1 

Butler- Volmer,  cathode  (Eq.  (13)) 

aa  =  1-4,  ac  =  0.6 

Exchange  current  density,  anode  (Eq. 

Xa  =  5.5x  1010  (Am-2),  a=  1, 

(18)) 

h-  1 ,  Eact.a  =  1 20  kJ  kmol- 1 

Exchange  current  density,  cathode  (Eq. 

Xc  =  7  x  109  (Am-2),  c  =  1, 

(19)) 

£act,c  =  120  kJ  kmol- 1 

Methane  reforming  (Eq.  (10)) 

&ref  =  8542  mol  bar-0  5  s- 1  m-2 

RT  (  1 

A  Vcon,c  —  “ I  “ q- 

neF  V-4 


iRTu 


exp 


2nQFDef  c  p 


AVcon,a  — 


RT 

nPF 


i - 


iRt  a 


In 


IFD^x^p 


\ 


iRt a 


\  1  +  2FDef,a*H20 P  ) 


(24) 


5.4.  Mass  balance 


Once  ideal  voltage  and  polarizations  are  known,  operating 
voltage  throughout  the  cell  can  be  evaluated  from  Eq.  (25) 

V  =  E  —  AFaet  -  AVohm  -  AVcon  (25) 

Given  the  operating  voltage  and  temperature  field  of  the  cell,  the 
latter  being  part  of  the  solution  to  the  thermal  model,  there  is 
only  one  value  of  current  density  satisfying  Eq.  (25)  at  each  slice. 
It  must  be  said  that  although  some  models  use  average  current 
density  along  the  cell  to  establish  operating  conditions  in  order 
to  diminish  computational  time,  this  do  not  correspond  to  real 
operation  of  a  fuel  cell.  Instead,  a  single  voltage  is  set  at  desired 
value  and  current  density  adjusts  to  satisfy  Eq.  (25)  according 
to  existing  temperatures  and  concentrations  inside  the  cell.  The 
model  presented  here  is  based  on  real  operating  conditions  so 
voltage  will  be  an  input  for  any  simulation. 

A  continuity  equation  must  be  applied  at  each  slice  after  solv¬ 
ing  Eq.  (25).  To  do  so,  variations  of  species  molar  flow  due 
to  methane  reforming,  water  gas  shift  and/or  hydrogen  oxida¬ 
tion  must  be  calculated.  Both  reforming  and  shift  are  evaluated 
according  to  Sections  2  and  4  with  values  shown  in  Table  4. 
Conversion  rate  of  hydrogen,  oxygen  and  water  due  to  electro¬ 
chemical  oxidation  is  calculated  thanks  to  Faraday’s  law 

iS 

Afiu2  = - -  =  —  A^h2o  =  2Aho2  (26) 

nQF 


stream  is 

/^outlet  ~  b iniet  A.hox[d  +  A/iref  +  (27) 

The  mass  balance  is  coupled  to  an  energy  balance  at  the  same 
control  volume,  which  comes  from  the  application  of  a  thermal 
model  to  heat  transfer  phenomena  inside  the  cell. 

6.  Thermal  model 

High  heat  flows  are  found  inside  an  SOFC  due  to  highly 
exothermic  and  endothermic  reactions  taking  place.  In  partic¬ 
ular,  the  heat  released  in  the  oxidation  of  hydrogen  is  mainly 
used  to  (i)  heat  the  air  flowing  through  the  air  supply  pipe  and 
(ii)  enhance  methane  reforming  at  the  anode. 

Heat  flows  must  be  controlled  due  to  mechanical  integrity 
issues  as  the  solid  structure  of  the  cell  is  formed  by  three  very  thin 
layers  whose  thicknesses  range  from  40  [xm  to  2  mm  and  with 
thermal  expansion  coefficients  between  10  and  13  x  l0-6  oC-l 
If  the  temperature  gradient  along  the  solid  tube  is  too  steep, 
thermal  expansion  mismatches  will  lead  to  poor  performance 
and,  eventually,  to  mechanical  failure.  For  these  reasons,  it  is 
necessary  to  have  a  highly  detailed  multidimensional  thermal 
model  in  order  to  simultaneously  search  for  high  efficiency  and 
reliability. 

The  following  heat  transfer  mechanisms  has  been  considered 
and  applied  to  control  volumes  in  Fig.  4: 

1 .  Axial  heat  conduction  inside  the  solid  structure  of  both  cell 
and  air  supply  pipe. 

2.  Radial  heat  convection  at  free  surfaces  including  inner  and 
outer  walls  of  cell  and  air  supply  pipe. 

3.  Radial  heat  radiation  between  cathode  inner  and  air  supply 
pipe  outer  walls. 

4.  Radial  heat  radiation  between  anodic  gas  and  anode  outer 
wall. 

6.1.  Conduction 


Finally,  the  expression  of  the  principle  of  conservation  of  mass  Solid  structures  are  considered  to  be  isothermal  in  the  radial 

to  be  applied  at  each  slice  and  for  each  component  of  each  direction  due  to  their  negligible  thickness.  Axially,  Fourier’s  law 
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is  used  to  determine  the  heat  flow  per  unit  area  perpendicular  to 
the  flow 


^cond  —  ^cd 


dT 

dx 


AT 

Ax 


(28) 


Estimating  thermal  conductivity  k  for  a  multilayer  structure  is  a 
difficult  task.  However,  a  first  guess  can  be  done  by  considering 
the  average  conductivity  as  proportional  to  the  conductivity  and 
thickness  of  each  layer  (Eq.  (29)) 


kc d,c4  T  ^cd,e4  T  ^cd,a4 


4  +  4  +  4 


(29) 


This  approach  is  not  rigorous  because  the  average  thermal  con¬ 
ductivity  should  be  closer  to  that  of  the  layer  acting  as  thermal 
barrier  or  insulator.  However,  an  analysis  of  ceramic  materi¬ 
als  typically  used  in  electrodes  and  electrolyte  manufacturing 
reveals  that  thermal  conductivity  ranges  from  2  to  4  W  K- 1  m_  1 , 
with  an  average  value  of  3  for  each  layer  [3] .  Thus,  the  error  made 
when  using  Eq.  (29)  to  evaluate  the  average  solid  structure  con¬ 
ductivity  is  negligible,  even  if  conductivity  is  different  from  one 
layer  to  another  (but  into  the  range  mentioned  before).  Finally, 
thermal  conductivity  of  air  supply  pipe  is  considered  to  be  a 
function  of  temperature  [1]. 


Fig.  9.  Cross-sectional  view  of  anodic  duct. 


6.2.  Convection 


An  accurate  study  of  convective  heat  transfer  is  absolutely 
necessary  if  maximum  temperature  of  electrodes  and  electrolyte 
are  not  to  be  exceeded.  It  must  be  noted  that  mean  operating 
temperature  and  maximum  temperature  are  very  close,  say  900 
and  1 100  °C,  respectively.  Besides,  convective  heat  transfer  phe¬ 
nomena  at  thermal  entries  may  be  the  cause  of  steep  temperature 
gradients  and,  consequently,  thermal  stresses.Newton’s  law  of 
cooling  has  been  used  to  evaluate  the  heat  flow  per  unit  area 


^cond  —  ^cv(2sup 


(30) 


where  hcv  is  the  convective  heat  transfer  or  film  coefficient. 

The  procedure  to  obtain  the  film  coefficient  at  each  one  of  the 
surfaces  is  quite  similar.  Firstly,  physical  properties  of  gases  in 
contact  with  each  surface  are  calculated:  specific  heat,  viscosity 
and  thermal  conductivity.  According  to  the  work  by  Todd  and 
Young  [18]  the  specific  heat  of  a  mixture  of  gases  is  proportional 
to  the  specific  heat  and  concentration  of  its  species.  The  specific 
heat  of  each  component  of  the  mixture  depends  on  temperature 
as  a  fifth  order  polynomial  while  the  dependence  of  viscosity  and 
thermal  conductivity  on  temperature  and  composition  is  more 
complex.  In  particular,  Reichenberg’s  and  Wassiljeva’s  expres¬ 
sions  are  used  to  evaluate  composition  and  thermal  conductivity, 
respectively,  although  the  scarce  data  for  conductivity  at  high 
temperature  may  lead  to  not  negligible  errors. 

In  second  place,  Reynolds  number  is  determined  locally,  for 
each  stream  and  each  slice.  For  this  purpose,  it  is  very  important 
that  the  geometry  of  each  duct  be  properly  defined,  specially 
the  anodic  duct  whose  complex  geometry  is  shown  in  Fig.  9.  In 
this  case,  the  hydraulic  diameter  is  calculated  from  the  general 


expression 

cross  section 

Ai  =  4 — - — (31) 

wet  perimeter 

Next  step  after  calculating  Reynolds  number  is  determining  the 
regime  of  heat  transfer  inside  the  duct.  Three  cases  are  consid¬ 
ered: 

1.  Faminar  heat  transfer,  thermal  and  hydrodynamic  entry 
length. 

2.  Faminar  heat  transfer,  fully  developed  flow. 

3.  Turbulent  heat  transfer,  fully  developed  flow. 

Thermal  and  hydrodynamic  boundary  layers  are  considered 
to  develop  at  the  same  time  as  long  as  Prandtl  number  is  near 
to  unity,  above  0.7.  The  critical  Reynolds  number,  based  on 
duct  hydraulic  diameter,  is  taken  as  2300  [20]  and  the  thermal 
entry  length  for  laminar  flow  is  given  by  the  following  equation 
[20,21]: 

jcte  =  0.0515RePrDh  (32) 

Where  the  flow  is  turbulent,  the  thermal  entry  length  is  taken  as 
negligible  according  to  [21]. 

Average  Nusselt  numbers  are  usually  evaluated  from  Re ,  Pr 
and  duct  geometry  through  empirical  expressions.  Nevertheless, 
this  approach  is  not  used  in  the  present  model  as  it  is  not  accu¬ 
rate  for  a  local  heat  transfer  analysis.  Instead,  correlations  for 
local  Nusselt  number  Nu  are  used.  Values  of  Nu  in  the  ther¬ 
mal  entry  region  are  given  by  Rohsenhow  et  al.  [20]  while  a 
constant  value  of  3.657  is  assumed  for  fully  developed  laminar 
flow.  Finally,  Gnielinski  equation  for  Nu  is  used  where  flow  is 
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turbulent  because  of  its  great  advantages  [21].  Usually,  corre¬ 
lations  for  forced  convection,  internal  flow  are  not  applicable 
to  transition  flows,  i.e.  Re  ranging  from  critical  value  to  about 
10,000,  due  to  uncertainty  about  flow  behaviour  but,  in  spite  of 
this  and  according  to  Baehr  and  Stephan  [21],  Gnielinski  equa¬ 
tion  is  valid  for  Re  above  2300.  Besides,  this  correlation  is  more 
accurate  with  respect  to  experimental  data  than  others,  which 
are  much  simpler,  like  Dittus  Boelter’s. 

The  convective  heat  transfer  coefficient  can  be  easily  evalu¬ 
ated  from  Nu 

krA  crcjc  -A llA 

hcv  =  cc'-as  (33) 


Table  5 

Coefficients  for  Hadvig’s  equation  [27] 

I  ki  bu  x  101  b2,i  x  104  bXi  x  107  b4,i  x  1011 


co2 


1 

0.3966 

0.4334 

2.620 

-1.560 

2.565 

2 

15.64 

-0.4814 

2.822 

-1.794 

3.274 

3 

394.3 

0.5492 

0.1087 

-0.3500 

0.9123 

h2o 

1 

0.4098 

5.977 

-5.119 

3.042 

-5.564 

2 

6.325 

0.5677 

3.333 

-1.967 

2.718 

3 

120.5 

1.800 

-2.334 

1.008 

-1.454 

C=f(T) 

=  b\  +  b2T+bT,T2 

+  b4T3 

(T  in  K) 

6.3.  Radiation 

The  analysis  of  radiative  heat  transfer  is  compulsory  if  an 
accurate  fuel  cell  model  is  to  be  developed,  the  main  reasons 
for  such  a  statement  being  a  high  view  factor  between  injector 
tube  and  cathode  and  the  high  temperatures  of  these  elements. 
Radiation  inside  the  cell  has  two  different  locations:  annular  duct 
between  cathode  and  air  injector  tube  and  anodic  duct. 

6.3.1.  Radiation  exchange  between  grey  walls 

In  the  annular  duct  radiative  exchange  between  grey  surfaces 
takes  place,  from  the  inner  face  of  the  cathode  to  the  outer  side  of 
the  injector  tube.  This  phenomenon  inside  fuel  cells  has  already 
been  included  in  some  of  the  models  developed  by  other  authors 
before  [9,13,22-24]  and  no  innovations  are  introduced  here. 
Grey  walls  satisfying  Kirchhoff ’s  laws  are  considered.  View  fac¬ 
tors  are  easily  evaluated  and  emissivity  is  taken  as  0.9  for  both 
surfaces.  The  total  heat  flow  exchanged  by  radiation  between 
cathode  and  air  supply  pipe  at  each  slice  is 

~  ^cathode  ^injector)  ~  A 

2 rad  —  “  ~  7  “  x- ^T^cathode  Ax  (34) 

1  |  r>>cathode  j  1  _  Y  j 

^cathode  ^injector  \  ^injector  J 

Radiative  exchange  between  walls  is  not  considered  in  the  anodic 
duct  (Fig.  9)  because,  due  to  symmetric  boundary  conditions,  all 
walls  are  considered  to  be  at  the  same  temperature.  Radiation 
heat  transfer  does  exist  but  the  net  heat  flux  at  each  wall  is  zero. 

6.3.2.  Grey  gas  radiation 

Radiative  heat  transfer  between  a  gas  volume  and  a  wall  is 
usually  negligible  except  if  the  optical  thickness  of  the  volume  is 
very  high  or  if  some  particular  gases  like  carbon  dioxide  and/or 
water  steam  are  present.  For  the  scope  of  this  work,  the  latter 
condition  is  satisfied  in  the  anode  as  the  fuel  gas  has  up  to  75% 
(v)  of  water  and  carbon  dioxide.  The  molar  fraction  of  these 
species  in  the  cathodic  gas  is  negligible  and  this  consideration 
does  not  apply  to  radiation  in  the  annular/cathodic  duct. 

Radiation  between  gases  and  walls  has  already  been  studied 
by  some  authors  before  [25,26]  but  these  works  present  some 
lack  of  accuracy  in  the  theoretical  treatment  of  the  gas  properties. 
VanderSteen  and  Pharoah  [25]  used  a  CFD  tool  to  analyze  the 
influence  of  both  kinds  of  radiation  on  the  simulation.  Results 
showed  that  radiation  between  walls  must  always  be  consid¬ 
ered.  With  respect  to  participating  media,  they  concluded  that 


their  influence  might  not  be  negligible  under  certain  operating 
conditions.  Samuelsen  [26]  included  the  effect  of  grey  gases  in 
all  ducts. 

Both  of  these  works  use  simplified  methods  for  estimating 
properties  of  gases  but  none  of  them  with  sufficient  accuracy. 
VanderSteen  and  Pharoah  studied  four  cases  ranging  from  no 
heat  absorption  from  gases  to  gases  with  high  absorption  coef¬ 
ficients  but  did  not  study  the  properties  of  gases  as  functions  of 
composition,  temperature  and  geometry.  Samuelsen  considered 
the  dependence  of  properties  on  composition  and  temperature 
using  the  work  by  Hottel  and  Sarofim  [20,21]  and  assumed  that 
gases  satisfy  Kirchhoff ’s  laws  in  order  to  simplify  equations. 

In  this  work,  properties  of  gases  are  determined  more  rigor¬ 
ously.  Emissivity  and  absorptivity  of  the  anodic  gas  depend  on 
its  composition,  temperature  and  on  the  geometry  of  the  duct 
through  the  mean  beam  length,  being  Hottel ’s  graphical  method 
the  easiest  way  to  calculate  them.  Nevertheless,  Hottel’s  method 
is  not  accurate  enough  and  Hadvig’s  analytical  method  for  cal¬ 
culating  the  emissivity  of  a  gas  mixture  [27]  is  used 


£gas  =  £C02  +  £H20  -  £C02£H20  (35) 

The  emissivity  of  each  gas  is  evaluated  from  the  following  equa¬ 
tion: 

3 

Si  =  yfCid  -  exp  (-kiPiL))  (36) 

i= 1 


where  pi  is  the  gas  partial  pressure  and  L  the  mean  beam  length. 
Coefficients  for  Eq.  (36)  can  be  found  in  Table  5. 

The  following  equation  given  in  [28]  is  used  to  evaluate  the 
absorptivity  of  the  gas  considering  that  gases  do  not  obey  Kirch¬ 
hoff’ s  laws 


^gas  — 


(  ^gas  A 
£co2  I  - -  ) 

V  1  surface  J 


0.65 


+  £H20 


f  ^gas  \ 
V  ^surface  J 


£co2£h2o 

(37) 


Using  these  radiative  properties  of  the  anodic  gas,  the  heat  flow 
exchanged  between  gas  and  surface  per  unit  area  can  be  calcu¬ 
lated 


9fad  — 


^surface0'" 


1  (1  0^gas)(l  £  surface) 


(^gas  ^gas  ^gas  ^surface  ) 


(38) 


where  all  the  walls  of  the  anodic  duct  are  supposed  to  be  at  the 
same  temperature,  as  said  before. 
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6.4.  Energy  balance 

Solving  energy  balance  equations  at  each  volume  of  each  slice 
gives  the  temperature  distribution  inside  the  cell.  The  following 
phenomena  must  be  taken  into  account: 

1.  Molar  mass  variations  due  to  gas  flow  inside  ducts. 

2.  Heat  release  and/or  consumption  due  to  reforming,  shift 
and/or  oxidation  reactions. 

3.  Heat  transfer  through  conduction,  convection  and  radiation. 

4.  Heating  due  to  Joule  effect  at  electric  current  conductors,  i.e. 
cathode,  anode  and  electrolyte. 

The  solution  to  the  electrochemical  model  is  used  as  an  input  to 
the  thermal  model  in  order  to  obtain  distributions  of  temperature, 
heat  flow  and  heat  transfer  coefficients  as  shown  in  Fig.  6.  More 
information  can  be  found  in  [19]. 

7.  Validation 

The  validation  of  the  model  is  a  two  step  process.  First,  a 
validation  of  the  global  performance  of  the  cell  is  done  by  com¬ 
paring  V-I  curves  from  tests  with  those  obtained  by  the  model. 
This  comparison  (Fig.  10)  shows  good  agreement  between  data 
published  by  Siemens-Westinghouse  from  tests  on  atmospheric 
cells  [29]  and  results  of  the  model  presented.  The  model  is  there¬ 
fore  capable  of  predicting  the  performance  of  the  tubular  fuel 
cell  with  accuracy. 

Besides  the  validation  made  in  Fig.  10,  a  numerical  evaluation 
of  the  model  has  been  done.  In  particular,  results  from  the  model 
at  atmospheric  and  overpressure  operation  have  been  compared 
with  those  in  reference  [1]  which  correspond  to  available  data 
from  tests  published  elsewhere.  This  comparison  is  shown  in 
Table  6  along  with  the  error  made  for  each  parameter. 

According  to  errors  shown  (Table  6),  which  are  smaller  than 
5%  for  the  most  important  data,  the  model  developed  is  accurate. 
In  addition,  it  is  remarkable  that  errors  decrease  considerably 
when  pressure  rises  because  this  means  that  the  model  can  be 
applied  to  fuel  cell-gas  turbine  hybrid  systems  with  minimal 
deviations  from  real  performance.  This  fact,  along  with  curve 


Fig.  10.  Comparison  with  data  from  Siemens-Westinghouse  [29]. 


Table  6 

Validation  (p  =  1.05/3.5  bar) 


Test 

Simulation 

Error  (%) 

Voltage  (V) 

0.69/0.64 

0.69/0.64 

— 

Current  density  (Am-2) 

1800/3000 

1905/2926 

5. 8/2.5 

Power  (W) 

104.8/157 

109.6/157 

4.6/0.5 

Fuel  utilization  (%) 

69/69 

71.3/68.2 

3.3/1. 1 

Air  utilization  (%) 

17.8/23.8 

18.6/23.4 

4.5/1. 7 

fitting  in  Fig.  10,  suggests  that  the  accuracy  of  the  model  pre¬ 
sented  above  is  satisfactory. 

A  direct  experimental  validation  of  the  thermal  model  is  not 
possible  as  long  as  there  are  no  data  available  for  the  temperature 
field  inside  the  cell.  Only  tests  conducted  by  Hirano  et  al.  [30] 
show  temperature  measurements  at  three  relevant  points  of  a 
tubular  cell:  initial,  middle  and  end  points.  However,  results  are 
not  comparable  since  the  cell  tested  by  Hirano  et  al.  is  a  30  cm 
long  cell  fed  mainly  with  hydrogen,  just  1%  methane.  Instead, 
a  comparison  with  a  previous  model  by  Campanari  and  Iora 
[1]  is  presented  in  order  to  evaluate  the  accuracy  of  the  model 
developed.  The  reference  model  is  similar  to  that  presented  here 
with  two  main  differences:  the  reforming  process  is  assumed 
to  be  kinetically  controlled  and  gas  radiation  is  not  considered. 
Figs.  11  and  12  show  temperature  and  anodic  gas  composition 
distributions  obtained  by  both  models;  similarity  is  observed  but 


Fig.  12.  Comparison  of  anode  gas  composition,  data  from  ref.  [1]. 
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Table  7 


Effect  of  anodic  gas  radiation 


Gas  radiation 

No  gas  radiation 

Voltage  (V) 

0.6 

0.6 

Current  density  (Am-2) 

3002 

3003 

Fuel  utilization  (%) 

80 

80 

Air  utilization  (%) 

19 

19 

Power  (W) 

150.2 

150.2 

Efficiency  (%) 

45.8 

45.8 

Mean  temperature  (°C) 

896 

897 

Peak  temperature  (°C) 

1003 

1003 

Case  1 :  0.6  V,  Uf  =  80%,  p-  3  bar. 

Table  8 

Effect  of  anodic  gas  radiation 

Gas  radiation 

No  gas  radiation 

Voltage  (V) 

0.2 

0.2 

Current  density  (Am-2) 

8031 

8097 

Fuel  utilization  (%) 

60.2 

60.7 

Air  utilization  (%) 

8.8 

8.9 

Power  (W) 

134 

135.1 

Efficiency  (%) 

11.5 

11.6 

Mean  temperature  (°C) 

599 

903 

Peak  temperature  (°C) 

1014 

1018 

Case  2:  0.2  V,  Uf  =  60%,  p  =  3  bar. 


not  coincidence.  The  differences  are  thought  to  be  caused  by  the 
improvements  made  as  they  are  located  in  the  first  part  of  the 
cell  where  reforming  is  at  equilibrium  (Figs.  7  and  8). 


8.  Effect  of  anodic  gas  radiation 

Two  main  subjects  have  been  introduced  in  this  model  with 
respect  to  previous  works:  a  mixed  method  to  evaluate  the  rate 
of  reforming  and  the  anodic  gas  radiation.  The  former  has  been 
studied  in  Section  4  and  its  influence  has  been  evaluated  both 
quantitatively  and  qualitatively.  The  effect  of  including  radiation 
between  anodic  gas  and  anode  outer  surface  is  now  analyzed. 

Two  cases  are  considered,  high  fuel  utilization  and  low  cur¬ 
rent  density  and  vice  versa;  pressure  is  set  at  3  bar  for  both  cases. 
First,  operation  at  0.6  V  and  80%  fuel  utilization  is  analyzed  for 
anodic  gas  taking  part  in  the  radiative  heat  transfer  and  not.  The 
results  of  both  simulations  are  shown  in  Table  7  and  it  is  con¬ 
cluded  that  for  low  current  density  the  effect  of  grey  gas  radiation 
is  negligible. 

Table  8  shows  results  for  the  second  case,  0.2  V  and  60% 
fuel  utilization.  Now,  results  do  not  coincide  when  gas  radia¬ 
tion  is  considered  or  not  as  including  the  effect  of  gas  radiation 
leads  to  poorer  performance  by  around  1%.  The  reason  for  this 
behaviour  is  the  increased  fuel  mass  flow.  When  current  density 
is  high,  fuel  mass  flow  increases,  the  anodic  gas  has  more  heat 
absorption  capacity  and  its  effect  as  heat  sink  is  higher,  result¬ 
ing  in  a  drop  in  cell  temperature.  This  decrease  in  temperature 
affects  to  the  electrochemical  performance  of  the  cell  negatively 
as  long  as  overpotential  increase,  mainly  ohmic  and  activation 
overpotentials. 


The  conclusion  to  the  comparison  shown  is  that  anodic  gas 
radiation  has  a  negligible  effect  on  performance  under  normal 
operating  conditions.  Fuel  cells  are  usually  operated  at  high 
voltage  and  low  current  density  seeking  for  high  efficiency 
although  it  may  be  interesting  to  increase  current  density  if  high 
power  density  is  sought,  despite  the  loss  in  efficiency.  However, 
if  fuel  flow  increases,  the  effect  of  gas  radiation  may  not  be 
negligible  and  should  be  considered.  This  situation  is  usually 
related  to  high  pressure,  low  voltage  and  low  fuel  utilization 
(Table  8). 

9.  Sensitivity  analysis  to  fuel  utilization  factor 

The  model  developed  allows  for  a  full  sensitivity  analysis  of 
all  the  parameters  determining  the  operating  conditions  of  the 
cell.  In  particular,  the  effect  of  fuel  utilization  shown  below  is 
part  of  a  complete  analysis  that  can  be  found  in  [19].  The  fuel 
utilization  factor  of  a  fuel  cell  accounts  for  the  fraction  of  fuel 
which  feeds  the  anode  and  is  used  to  produce  electrical  work 
(Eq.  (39)).  It  must  not  be  taken  as  a  sort  of  efficiency 

Uf  =  %2lused -  (39) 

W  H2 ,  equivalent  I  inlet 

The  influence  of  fuel  utilization  on  cell  performance  is  signif¬ 
icant  as  depicted  in  Fig.  13,  which  shows  the  voltage-current 
density  curve  for  two  values  of  this  parameter.  It  is  important  to 
note  that  the  effect  of  increasing  fuel  utilization  is  not  the  same 
at  high  or  low  current  densities. 

An  analysis  of  fuel  cell  internal  performance  has  been  done 
in  order  to  understand  Fig.  13.  The  Nernst  potential  depends 
on  temperature  and  composition  as  shown  above  (Eqs.  (11)  and 
(12)).  Figs.  14-16  depict  the  distribution  of  these  parameters 
along  the  cell,  temperature  being  drawn  only  for  the  solid  struc¬ 
ture. 

Nernst  potential  increases  with  hydrogen  and  oxygen  con¬ 
centrations  and  decreases  with  water  concentration  (Eq.  (11)). 
Thus,  according  to  Figs.  14  and  15,  the  potential  at  90%  fuel  uti¬ 
lization  should  be  higher  for  the  first  50  cm  of  the  cell  despite  the 
effect  of  oxygen  concentration  as  this  is  diminished  by  a  power 
to  0.5.  For  the  second  half  of  the  cell,  the  potential  should  be 


Fig.  13.  Effect  of  fuel  utilization  on  voltage-current  density  curve. 
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Fig.  14.  Influence  of  fuel  utilization  on  anodic  gas  composition,  V=  0.6  V. 


Fig.  15.  Influence  of  fuel  utilization  on  oxygen  concentration  at  cathode, 
V=0.6V. 

higher  for  60%  fuel  utilization  as  concluded  from  compositions 
(Figs.  14  and  15). 

The  effect  of  fuel  utilization  on  temperature  is  depicted  in 
Fig.  1 6  where  it  can  be  seen  that  decreasing  this  parameter  causes 
initial  subcooling  and  final  overheating  with  respect  to  normal 
values  around  80%.  When  the  fuel  utilization  factor  is  low,  keep¬ 
ing  voltage  constant,  there  is  a  high  fraction  of  fuel,  which  is  not 
used  for  producing  power.  This  fraction  of  fuel  is  not  oxidized 


Fig.  16.  Influence  of  fuel  utilization  on  temperature  of  solid  structure,  V  =  0.6  V. 


Fig.  17.  Influence  of  fuel  utilization  on  Nernst  potential,  V=0.6  V. 


but  it  is  reformed,  originating  a  significant  decrease  in  temper¬ 
ature  in  the  first  part  of  the  cell.  The  increase  in  temperature  in 
the  second  half  of  the  cell  is  caused  by  a  lack  of  cooling  air,  con¬ 
sidering  that  the  average  operating  temperature  and  boundary 
conditions  for  natural  gas  and  air  streams  feeding  the  cell  are 
the  same  for  all  cases.  The  Nernst  potential  is  directly  affected 
by  these  changes  in  temperature  (Eq.  (12)),  decreasing  where 
temperature  is  higher. 

It  can  be  concluded  that  the  effect  of  temperature  overcomes 
that  of  composition  (Fig.  17).  Thus,  ideal  voltage  (Fig.  17)  is 
higher  where  temperature  is  higher  (Fig.  16)  despite  a  delay 
in  the  cross  point  of  voltage  lines  with  respect  to  temperature 
lines  which  take  place  at  50  and  80  cm  from  the  entrance  to 
the  cell,  respectively.  This  delay  is  caused  by  the  effect  of  gas 
composition,  whose  influence  opposes  that  of  temperature. 

10.  Conclusions 

The  potential  and  reliability  of  the  model  developed  has  been 
shown.  This  computational  tool  includes  a  very  detailed  electro¬ 
chemical  model  of  performance,  which  is  the  result  of  a  com¬ 
plete  analysis  of  previous  works  and  includes  new  approaches 
to  critical  tasks  such  as  internal  reforming  of  natural  gas.  In 
this  case,  reaction  equilibrium  and  kinetics  have  been  consid¬ 
ered. 

The  simulation  of  heat  transfer  phenomena  inside  the  cell  has 
been  done  in  a  very  careful  way  so  that  no  information  about 
local  heat  flows  and/or  temperatures  were  lost  or  neglected.  For 
this  purpose,  radiation  between  anodic  gas  and  surface  has  been 
considered. 

Particular  conclusions  are  listed  below: 

1.  The  effect  of  anodic  gas  radiation  can  be  neglected  under 
normal  operating  conditions.  However,  results  should  be  low¬ 
ered  by  1  %  if  fuel  flow  increases  considerably,  this  situation 
being  related  to  high  current  density  and  pressure  and  low 
fuel  utilization  factor. 

2.  The  behaviour  of  the  reforming  reaction  is  not  uniform  along 
the  cell.  Care  should  be  taken  when  considering  whether 
equilibrium  is  reached  or  not,  specially  when  methane  con- 
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centration  drops  down  and  temperature  increases.  These  con¬ 
ditions  are  usually  found  at  the  second  half  of  the  cell. 

3.  The  model  presented  applies  to  cells  in  the  inner  part  of  the 
stack  but  a  fully  three-dimensional  study  must  be  done  in 
order  to  evaluate  the  influence  of  peripheral  cells  on  the  global 
performance  of  the  stack. 

It  can  be  concluded  that  the  tool  presented  can  help  to  develop 
tubular  solid  oxide  fuel  cell  technology,  improving  both  cell  per¬ 
formance  and  reliability.  Multidimensional  models  are  essential 
to  fuel  cell  researchers  and  manufacturers  as  long  as  installation 
costs  are  still  very  high  and  research  cannot  rely  on  tests  entirely. 
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